Devel_int WGCNA modules expression in NT D5, D7, and D10 samples

Render this report with ~/spinal_cord_paper/scripts/Gg_devel_scWGCNA_module_analysis_render.sh.

library(Seurat)

Attaching SeuratObject

library(WGCNA)

Loading required package: dynamicTreeCut
Loading required package: fastcluster

Attaching package: ‘fastcluster’

The following object is masked from ‘package:stats’:

    hclust



Attaching package: ‘WGCNA’

The following object is masked from ‘package:stats’:

    cor

library(tidyr)
library(ggplot2)
library(stringr)
library(patchwork)
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ readr   1.4.0      ✔ forcats 0.5.1 
✔ purrr   0.3.4      
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

library(cowplot)

Attaching package: ‘cowplot’

The following object is masked from ‘package:patchwork’:

    align_plots

library(pheatmap)
library(gridExtra)

Attaching package: ‘gridExtra’

The following object is masked from ‘package:dplyr’:

    combine

source("~/Neuraltube/scripts/heatmap4.R")

Load individual seurat and test WGCNA data

The individual data sets are the Day 5 (Gg_D05_ctrl_seurat_070323), Day 7 (Gg_D07_ctrl_seurat_070323), and Day 10 (Gg_ctrl_1_seurat_070323) chicken spinal cord sets. The test WGCNA data are the modules calculated on the integrated data set of all three stages.

se_path <- c("Gg_D05_ctrl_seurat_070323",
             "Gg_D07_ctrl_seurat_070323",
             "Gg_ctrl_1_seurat_070323")

clust_col <- read.csv("~/spinal_cord_paper/annotations/broad_cluster_marker_colors.csv") %>% 
  rename(broad = broad_cluster) %>% 
  select(-marker)

Order of the broad clusters for plotting purposes.

broad_order <- c("progenitors",
      "FP",
      "RP",
      "FP/RP",
      "neurons",
      "OPC",
      "MFOL",
      "pericytes",
      "microglia",
      "blood",
      "vasculature"
      )

my.ses <- list()
col_table <- list()
ord_levels <- list()

for (i in seq(se_path)) {
  # load the data sets
  my.se <- readRDS(paste0("~/spinal_cord_paper/data/", se_path[i], ".rds"))
  annot <- read.csv(list.files("~/spinal_cord_paper/annotations",
                               pattern = str_remove(se_path[i], "_seurat_\\d{6}"),
                               full.names = TRUE))
  
  if(length(table(annot$number)) != length(table(my.se$seurat_clusters))) {
     stop("Number of clusters must be identical!")
  }
  
  # rename for left join
  annot <- annot %>% 
    mutate(fine = paste(fine, number, sep = "_")) %>% 
    mutate(number = factor(number, levels = 1:nrow(annot))) %>% 
    rename(seurat_clusters = number) 
  
  # cluster order for vln plots
  ord_levels[[i]] <- annot$fine[order(match(annot$broad, broad_order))]
  
  # create index for color coding
  col_table[[i]] <- annot %>%
    left_join(clust_col, by = "broad") %>% 
    select(c("fine", "color"))
  
  # add cluster annotation to meta data
  my.se@meta.data <- my.se@meta.data %>% 
    rownames_to_column("rowname") %>% 
    left_join(annot, by = "seurat_clusters") %>% 
    mutate(fine = factor(fine, levels = annot$fine)) %>% 
    column_to_rownames("rowname")
  
  my.ses[[i]] <- my.se

}

names(my.ses) <- c("Gg_D05_ctrl", "Gg_D07_ctrl", "Gg_D10_ctrl")
names(col_table) <- c("Gg_D05_ctrl", "Gg_D07_ctrl", "Gg_D10_ctrl")
names(ord_levels) <- c("Gg_D05_ctrl", "Gg_D07_ctrl", "Gg_D10_ctrl")

rm(my.se, annot)

# The reference WGCNA data. We can have several, if we want to test many at the same time
WGCNA_data = list()
WGCNA_data[[1]] = readRDS("~/spinal_cord_paper/output/Gg_devel_int_scWGCNA_250723.rds")
my.wsub =list()
my.wsub[[1]]= c(1:22)

# the name of each sample, as they appear in my.files and in the metadata of the combined object
my.samplenames = c("D05", "D07", "D10")

# The subset of clusters in each of the corresponding samples
my.subset=list(c(1:24),
               c(1:26),
               c(1:22))

# This is just to add a little bit more sense to the modules, so that we don't get just a color. Corresponds to WGCNA_data
my.modulenames = list()
my.modulenames[[1]] = c(1:22)

#Subset the seurat objects if needed as defined above
for (i in 1:length(my.ses)) {
  my.ses[[i]] = subset(my.ses[[i]], idents = my.subset[[i]])
}

# Take only genes that are present in the all samples
my.dcs = list()

# We go trhough each WGCNA data
for(i in 1:length(WGCNA_data)) {
  #We start with complete modules
  my.dc = WGCNA_data[[i]]$dynamicCols
  
  #Remove the few genes that are not found in the other datasets
  my.dc = my.dc[which(names(my.dc) %in% rownames(my.ses[[1]]@assays$RNA@data))]
  my.dc = my.dc[which(names(my.dc) %in% rownames(my.ses[[2]]@assays$RNA@data))]
  my.dc = my.dc[which(names(my.dc) %in% rownames(my.ses[[3]]@assays$RNA@data))]
  
  my.dc = my.dc[which(my.dc %in% names(table(WGCNA_data[[i]]$dynamicCols))[my.wsub[[i]]])]
  
  my.dcs[[i]] = my.dc
  
}

Module gene correlation

Here, we do a correlation matrix / heatmap, to see which cell clusters group togheter. This helps us to make more detailed dotplots.
This part of the script can still be used to compare several WGCNA datasets in parallel.

# broad cluster color table
all_col <- do.call(rbind, col_table) %>% 
  rownames_to_column("sample") %>% 
  mutate(sample = substr(sample, 1, 3)) %>% 
  mutate(sample_celltype = paste(sample, fine, sep = "_")) %>% 
  select(c("color", "sample_celltype", "sample"))

all_col$sample

 [1] "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05" "D05"
[25] "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07" "D07"
[49] "D07" "D07" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10" "D10"

spearman correlation heatmap

annotations

# names and colors for the heatmap annotation
annot_name <- data.frame(
  "Celltypes" = all_col$sample_celltype,
  "Sample"    = all_col$sample,
  row.names = all_col$sample_celltype)
  

pheat_col_table <- do.call(rbind, col_table) %>% 
  rownames_to_column("sample") %>% 
  mutate(sample = substr(sample, 1,11)) %>% 
  mutate(fine = paste(sample, fine, sep = "_"))

# match color table with annotation
pheat_col_table <- pheat_col_table[match(annot_name$Celltypes, pheat_col_table$fine),]

annot_col <- list(
  Celltypes = pheat_col_table$color,
  Sample = c(D05 = "#A4A4A4",
             D07 = "#515151",
             D10 = "#000000")
  )

names(annot_col[[1]]) <- annot_name$Celltypes

heatmap of module pseudobulk average expression

Load the integrated data set

The integrated data set on which the WGCNA is calculated. We plot it, split by the original cell types from the three original samples.

Avg. module exp. by stage

tSNE DimPlots showing the average expression of each module by stage.

AE over time

We plot the average expression of each module in the three stages and the 5 broad cell type clusters present in all 3 stages.

mod_annot <- read.csv("~/spinal_cord_paper/annotations/Gg_devel_int_scWGCNA_module_annotation.csv") %>%
  dplyr::mutate(module = str_replace_all(module, "\\d{1,2}\\_", "AE"))

meta <- list()

for (i in seq(my.ses)) {
  meta[[i]] <- my.ses[[i]]@meta.data %>%
    tibble::rownames_to_column("cell_ID") %>%
    dplyr::mutate(cell_ID = paste0(cell_ID, "_", i)) %>%
    dplyr::select(c("cell_ID", "broad"))
}

meta <- do.call(rbind, meta)

# mean average expression by stage
mean_AE <- avg.mod.eigengenes %>%
  tibble::rownames_to_column("cell_ID") %>%
  dplyr::mutate(stage = stringr::str_sub(cell_ID, -1)) %>%
  dplyr::mutate(stage = factor(stage, levels = c(1:3), labels = c("D05", "D07", "D10"))) %>%
  dplyr::left_join(meta, by = "cell_ID") %>%
  tibble::column_to_rownames("cell_ID") %>%
  tidyr::unite("stage_cl", stage, broad, sep = "_") %>%
  dplyr::group_by(stage_cl) %>%
  dplyr::summarise_each(mean) %>%
  dplyr::ungroup() %>%
  gather(key="module", value = "AE", -stage_cl) %>%
  dplyr::left_join(mod_annot[, c(1,3)], by = "module") %>%
  tidyr::separate("stage_cl", c("stage", "broad"), sep = "_", remove = FALSE) 

Warning: `summarise_each_()` was deprecated in dplyr 0.7.0.
ℹ Please use `across()` instead.
ℹ The deprecated feature was likely used in the dplyr package.
  Please report the issue at <https://github.com/tidyverse/dplyr/issues>.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

# %>%
  # dplyr::filter(broad %in% c("progenitors", "neurons", "RP", "FP", "pericytes"))

labels_dotplot <- stringr::str_remove(modules_in_order, "^AE")
names(labels_dotplot) <- modules_in_order

plot_clusters <- c("progenitors", "neurons", "RP", "FP", "pericytes", "OPC", "MFOL", "microglia", "blood")

mean_mod <- ggplot(data = mean_AE,
  aes(
    x = stage,
    y = AE,
    color = factor(broad, levels = plot_clusters),
    group = broad,
    label = annotation
    )
  ) +
  geom_line() +
  geom_point() +
  scale_color_manual(values = clust_col$color[match(plot_clusters, clust_col$broad)]) +
  theme_bw() +
  # facet wrap with reordered factors
  facet_wrap(vars(factor(module, levels = unique(mean_AE$module)[module_order])),
             scales = "free_y",
             nrow = 4,
             ncol = 6) +
  labs(color = "broad") +
  ggtitle("Average module expression by stage")

plotly::ggplotly(mean_mod)

VlnPlots of avg. module exp. by stage and seurat cluster

colored by module

# reorder seurat clusters
for (i in seq(my.ses)) {
  my.ses[[i]]$seurat_clusters <- factor(
  my.ses[[i]]$seurat_clusters,
  levels = levels(my.ses[[i]]$seurat_clusters)[as.integer(str_extract(ord_levels[[i]], "\\d{1,2}$"))]
  )

}

vplots <- list()

for (i in seq(my.ses)) {
  
  mods <- colnames(my.ses[[i]][[]])[grep("^AE",colnames(my.ses[[i]][[]]))]
  
  vplots[[i]] <- VlnPlot(
        my.ses[[i]],
        features = mods[module_order],
        group.by = "seurat_clusters",
        stack = TRUE, flip = TRUE,
        cols = substring(mods, 3)[module_order]) +
    theme(legend.position = "none") +
    geom_hline(yintercept = 0, lty = "dashed")
  
}

vplots[[1]]
vplots[[2]]

vplots[[3]]

colored by cell type

clust_col <- read.csv("~/spinal_cord_paper/annotations/broad_cluster_marker_colors.csv")

vplots_id <- list()

for (i in seq(my.ses)) {
    
  mods <- colnames(my.ses[[i]][[]])[grep("^AE",colnames(my.ses[[i]][[]]))]
  
  vplots_id[[i]] <- VlnPlot(
    my.ses[[i]],
    features = mods[module_order],
    group.by = "seurat_clusters",
    stack = TRUE, flip = TRUE,
    fill.by = "ident",
    cols = col_table[[i]]$color[as.integer(str_extract(ord_levels[[i]], "\\d{1,2}$"))]) +
    theme(legend.position = "none") +
    geom_hline(yintercept = 0, lty = "dashed")
}

vplots_id[[1]]
vplots_id[[2]]

vplots_id[[3]]

MN modules

For the figures we specifically select the two MN modules and plot them as Vln plots.

avg. module exp. by stage

v1 <- VlnPlot(my.sec,
        features = mods[module_order], 
        group.by = "orig.ident")
v1

pdf("~/spinal_cord_paper/figures/Fig_1_AE_by_stage.pdf", height = 20, width = 20)
v1

module GO terms barplots

# Load the top 50 (by limma::topGO) Go term table list from scWGCNA_Gg_devel_int.Rmd
goterms <- readRDS("~/spinal_cord_paper/output/Gg_devel_int_module_GOTerms_250723.rds")
# Load the devel module annotations 
modules <- read.csv("~/spinal_cord_paper/annotations/Gg_devel_int_scWGCNA_module_annotation.csv")
# Set list names as module number and name (e.g. "1_black")
names(goterms) <- modules$module

plot_list <- list()

for (i in seq(goterms)) {
  # Select the darkred module (motor neuron/transmembrane transport)
  toplot <- goterms[[i]]
  terms <- as.character(toplot$Term)
  
  #Plot the barplot!
  plot_list[[i]] <- ggplot(toplot, aes(x=Term, y=-log10(P.DE), fill=Ont)) +
    geom_bar(stat="identity") +
    scale_fill_manual(values = str_remove(names(goterms)[i], "^\\d{1,2}_")) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1),
          axis.text.y = element_text(size = 5)) +
    labs(title=paste0("Top GOTerm of module ", names(goterms)[i])) +
    coord_flip() +
    geom_hline(yintercept = -log10(0.05), lty = "dashed") +
    theme_cowplot()
}

rm(toplot, terms)

names(plot_list) <- modules$module

# Export PDF
pdf("~/spinal_cord_paper/figures/GO_term_barplot_scWGCNA_Gg_devel_modules.pdf", width = 14, height = 7)
for (i in seq(goterms)) {
  grid.arrange(plot_list[[i]])
}
dev.off()

# Date and time of Rendering
Sys.time()

sessionInfo()